/*
 * numeric.cpp
 *
 *  Created on: Feb 23, 2012
 *      Author: falko
 */

#include "numeric.h"
#include <iostream>
#include <sstream>


namespace falko
{


hpf_t::hpf_t()
{
    mpfr_init2(val_, GMP_PREC);
    mpfr_set_zero(val_, 1);
}

hpf_t::hpf_t(double val)
{
    mpfr_init2(val_, GMP_PREC);
    mpfr_set_d(val_, val, rounding_type_);
}

hpf_t::hpf_t(hpf_t const &v)
{
    mpfr_init2(val_, GMP_PREC);
    mpfr_set(val_, v.val_, rounding_type_);

}

hpf_t::hpf_t(const char* v, uint base)
{
    mpfr_init2(val_, GMP_PREC);
    mpfr_set_str(val_, v, base, rounding_type_);
}

hpf_t hpf_t::operator=(hpf_t const &v)
{
	mpfr_set(val_, v.val_, rounding_type_);
    return *this;
}

hpf_t hpf_t::operator=(double v)
{
    mpfr_set_d(val_, v, rounding_type_);
    return *this;
}

hpf_t hpf_t::operator=(const char* v)
{
    mpfr_set_str(val_, v, 10, rounding_type_);
    return *this;
}

hpf_t &hpf_t::operator/=(hpf_t const &v)
{
    try
    {
        mpfr_div(val_, val_, v.val_, rounding_type_);
    }
    catch (...)
    {
        std::cerr << "GMP: divide by zero" << std::endl;
        mpfr_set_d(val_, 0.0, rounding_type_);
    }
    return *this;
}

hpf_t hpf_t::operator-() const
{
    hpf_t retVal;
    mpfr_neg(retVal.val_, retVal.val_, rounding_type_);
    return hpf_t(retVal);
}

hpf_t hpf_t::operator+(hpf_t const &v) const
{
    hpf_t retVal;
    mpfr_add(retVal.val_, val_, v.val_, rounding_type_);
    return hpf_t(retVal);
}

hpf_t hpf_t::operator-(hpf_t const &v) const
{
    hpf_t retVal;
    mpfr_sub(retVal.val_, val_, v.val_, rounding_type_);
    return hpf_t(retVal);
}

hpf_t hpf_t::operator*(hpf_t const &v) const
{
    hpf_t retVal;
    mpfr_mul(retVal.val_, val_, v.val_, rounding_type_);
    return hpf_t(retVal);
}

hpf_t hpf_t::operator*(double v) const
{
    hpf_t retVal;
    mpfr_mul(retVal.val_, val_, hpf_t(v).val_, rounding_type_);
    return hpf_t(retVal);
}

hpf_t hpf_t::operator/(hpf_t const &v) const
{
    hpf_t retVal = *this;
    retVal /= v;
    return hpf_t(retVal);
}

hpf_t hpf_t::operator/(uint v) const
{

    if (v == 0)
        return *this;
    hpf_t retVal;
    mpfr_div_ui(retVal.val_, val_, v, rounding_type_);
    return hpf_t(retVal);
}

hpf_t hpf_t::pow(uint v) const
{
    hpf_t retVal;
    mpfr_pow_ui(retVal.val_, val_, v, rounding_type_);
    return hpf_t(retVal);
}

hpf_t hpf_t::sqrt() const
{
    hpf_t retVal;
    mpfr_sqrt(retVal.val_, val_, rounding_type_);
    return hpf_t(retVal);
}

hpf_t hpf_t::abs() const
{
    hpf_t retVal;
    mpfr_abs(retVal.val_, val_, rounding_type_);
    return hpf_t(retVal);
}

bool hpf_t::operator==(hpf_t const &v) const
{
    return (mpfr_cmp(val_, v.val_) == 0);
}

bool hpf_t::operator!=(hpf_t const &v) const
{
    return !(mpfr_cmp(val_, v.val_) == 0);
}

bool hpf_t::operator<(hpf_t const &v) const
{
    return (mpfr_cmp(val_, v.val_) < 0);
}

bool hpf_t::operator<=(hpf_t const &v) const
{
    return (mpfr_cmp(val_, v.val_) < 0) ||
           (mpfr_cmp(val_, v.val_) == 0);
}

bool hpf_t::operator>(hpf_t const &v) const
{
    return (mpfr_cmp(val_, v.val_) > 0);
}

bool hpf_t::operator>=(hpf_t const &v) const
{
    return (mpfr_cmp(val_, v.val_) > 0) ||
           (mpfr_cmp(val_, v.val_) == 0);
}

double hpf_t::toDouble() const
{
    return mpfr_get_d(val_, rounding_type_);
}


void hpf_sin(hpf_t &rop, hpf_t const &op)
{
	mpfr_sin(rop.val_, op.val_, rop.rounding_type_);
}
void hpf_cos(hpf_t &rop, hpf_t const &op)
{
	mpfr_cos(rop.val_, op.val_, rop.rounding_type_);
}



hpv3d_t::hpv3d_t(hpf_t const &x, hpf_t const &y, hpf_t const &z)
{
    x_ = x;
    y_ = y;
    z_ = z;
}

hpv3d_t::hpv3d_t(hpv3d_t const &v)
{
    x_ = v.x_;
    y_ = v.y_;
    z_ = v.z_;
}

hpv3d_t hpv3d_t::operator= (hpv3d_t const &v)
{
    x_ = v.x_;
    y_ = v.y_;
    z_ = v.z_;
    return *this;
}

hpv3d_t  &hpv3d_t::operator+=(hpv3d_t const &v)
{
    x_ += v.x_;
    y_ += v.y_;
    z_ += v.z_;
    return *this;
}

hpv3d_t &hpv3d_t::operator-=(hpv3d_t const &v)
{
    x_ -= v.x_;
    y_ -= v.y_;
    z_ -= v.z_;
    return *this;
}

hpv3d_t &hpv3d_t::operator*=(hpf_t const &v)
{
    x_ *= v;
    y_ *= v;
    z_ *= v;
    return *this;
}

hpv3d_t &hpv3d_t::operator%=(hpv3d_t const &v)
{
    hpf_t x = (y_ * v.z_) - (z_ * v.y_);
    hpf_t y = (z_ * v.x_) - (x_ * v.z_);
    hpf_t z = (x_ * v.y_) - (y_ * v.x_);
    x_ = x;
    y_ = y;
    z_ = z;
    return *this;
}

hpv3d_t &hpv3d_t::operator/=(hpf_t const &v)
{
    x_ /= v;
    y_ /= v;
    z_ /= v;
    return *this;
}

hpv3d_t hpv3d_t::operator-() const
{
    return hpv3d_t(-x_, -y_, -z_);
}

hpv3d_t hpv3d_t::operator+(hpv3d_t const &v) const
{
    hpv3d_t retVal(*this);
    retVal += v;
    return hpv3d_t(retVal);
}

hpv3d_t hpv3d_t::operator-(hpv3d_t const &v) const
{
    hpv3d_t retVal(*this);
    retVal -= v;
    return hpv3d_t(retVal);
}

hpv3d_t hpv3d_t::operator%(hpv3d_t const &v) const
{
    hpv3d_t retVal(*this);
    retVal %= v;
    return hpv3d_t(retVal);
}

hpv3d_t hpv3d_t::operator*(hpf_t const &v) const
{
    hpv3d_t retVal(*this);
    retVal *= v;
    return hpv3d_t(retVal);
}

hpf_t hpv3d_t::operator*(hpv3d_t const &v) const
{
    return hpf_t(x_ * v.x_ + y_ * v.y_ + z_ * v.z_);
}

hpv3d_t hpv3d_t::operator/(hpf_t const &v) const
{
    hpv3d_t retVal(*this);
    retVal /= v;
    return hpv3d_t(retVal);
}

bool hpv3d_t::operator==(hpv3d_t const &v) const
{
    return x_ == v.x_ && y_ == v.y_ && z_ == v.z_;
}

bool hpv3d_t::operator!=(hpv3d_t const &v) const
{
    return x_ != v.x_ || y_ != v.y_ || z_ != v.z_;
}

void hpv3d_t::normalize()
{
    *this /= magnitude();
}

hpf_t hpv3d_t::magnitude() const
{
    return (x_.pow(2) + y_.pow(2) + z_.pow(2)).sqrt();
}

void HpFileWriter::writeHpf(hpf_t const &v)
{
//	std::cout << v.toDouble() << std::endl;
    mpfr_out_str(out_file_, 10, prec_, v.val_, MPFR_RNDNA);
//    fprintf(out_file_, "%lf", v.toDouble());
}

void HpFileWriter::writeHpv3d(hpv3d_t const &v)
{
    writeHpf(v.x_);
    writeStr("\t");
    writeHpf(v.y_);
    writeStr("\t");
    writeHpf(v.z_);
}

void HpFileWriter::writeStr(char const *v)
{
    fprintf(out_file_, "%s", v);
}

HpFileWriter& operator<<(HpFileWriter& sout, hpf_t const &v)
{
    sout.writeHpf(v);
    return sout;
}
HpFileWriter& operator<<(HpFileWriter& sout, hpv3d_t const &v)
{
    sout.writeHpv3d(v);
    return sout;
}
HpFileWriter& operator<<(HpFileWriter& sout, char const *v)
{
    sout.writeStr(v);
    return sout;
}

} // End Falko
